Theory of the Magnetic phase diagram of LiHoF 4 
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The properties of LiHoF4 are believed to be well described by a long-range dipolar Ising model. 
We go beyond mean-field theory and calculate the phase diagram of the Ising model in a transverse 
field using a quantum Monte Carlo method. The relevant Ising degrees of freedom are obtained 
using a non-perturbative projection onto the low-lying crystal field eigenstates. We explicitly take 
the domain structure into account, and the strength of the near-neighbor exchange interaction 
is obtained as a fitting parameter. The on-site hyperfine interaction is approximately taken into 
account through a renormalization of the transverse applied magnetic field. Finally, we propose a 
spectroscopy experiment to precisely measure the most important parameter controlling the location 
of the phase boundary. 

PACS numbers: 75.10. Jm,75.40.Mg,75.70.Ak,75.30.Gw 



I. INTRODUCTION 

In the last decade, the rare-earth compound LiHoF4 
has been found to display an array of interesting mag- 
netic phenomena. At high temperatures LiHoF4 is a 
paramagnet, but there is a second-order transition to 
a ferromagnetic state at 1.53 Ki This Ising magnetic 
transition is driven by the weak magnetic dipole inter- 
action, with a strength of order IK, and not the more 
usual Coulomb exchange interaction. The critical tem- 
perature can be lowered by application of a magnetic 
field transverse to the easy-axis direction of ferromag- 
netic ordering. The magnetic field introduces quantum 
fluctuations of the spins and beyond a critical value of 
~ 4.9 Tesla, destroys long-range ordering even at zero 
temperature. LiHoF4 thus represents a model magnet for 
studying quantum phase transitions^ Since the rate of 
quantum tunnelling between different spin configurations 
can be carefully controlled with the transverse magnetic 
field, this material constitutes a good testing ground for 
the efficiency of quantum annealing^ By substituting the 
magnetic Ho 3+ ions with non-magnetic Y 3+ ions, disor- 
der can be introduced, and spin-glass behavior has been 
observed when the magnetic ions are sufficiently dilute. 4 
On further dilution the range of dynamic time scales dis- 
plays a remarkable narrowing in what has been called the 
"antiglass" phased 

The magnetic properties of LiHoF4 originate in the 
Ho 3+ ions. The ground state of the Ho 3+ ion in the 
crystal field is a doublet, and the first excited state is 
~ 11K above the ground state. The crystal field states of 
the Ho 3+ ion are such that there are no matrix elements 
of the transverse angular momentum ( J Xl J y ) within the 
ground state doublet. Hence the transverse susceptibility 
vanishes (to lowest order in the the applied field) giving 
rise to strong Ising anisotropy. Therefore LiHoF4 in a 



transverse magnetic field, and at temperatures lower than 
~ 11K, is believed to be a very good realization of a 
dipolar Ising model 
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where J is the coupling constant, nj the interspin dis- 
tance and Zij the component of the interspin distance 
along the Ising axis. The effective transverse field pa- 
rameter h x is a measure of the (higher order) mixing 
effects introduced by the physical transverse magnetic 
field. The summation is done over all Ho 3+ ions, which 
sit on a body-centered tetragonal lattice with four Ho 3+ 
ions per unit cell^ 

The goal of the present study is to determine the quan- 
titative phase diagram of pure L1H0F4 from (quasi-) first 
principles, starting from a crystal-field hamiltonian that 
has been fit to spectroscopic data. Bitko et ali£ suc- 
cessfully fit their phase diagram data using a mean-field 
theory with two free parameters, a transverse suscepti- 
bility g± w 0.74 to replace the Lande g factor = 1.25 
which thus rescales the transverse field, and an effective 
dipole coupling strength Jo which rescales the tempera- 
ture. This calculation also did not take into account the 
domain structure of the ferromagnetic state. In another 
calculation, R0nnow et al£ have recently used an RPA 
method to find the collective mode softening seen in their 
neutron scattering measurements and obtain an estimate 
of the phase diagram. 

The advancements reported in the present study over 
earlier attempts are twofold. We develop an effective low- 
energy hamiltonian (derived non-perturbatively by pro- 
jection from the full 17 state crystal field hamiltonian) 
that acts on spin-i Ising degrees of freedom, and we go 
beyond simple mean-field theory by using extensive quan- 
tum Monte Carlo simulations. An earlier classical Monte 
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Carlo study found a critical temperature of 1.89 K by ex- 
trapolation from rather limited system sizesi However, 
they did not take into account a near-neighbor exchange 
interaction among the Ho 3+ ions and the structure of 
domains observed in L1H0F4. Yet another Monte Carlo 
study finds a transition temperature of 1.51 K, but only 
by adjusting the strength of the dipolar interaction to 
reproduce the experimentally determined ground-state 
energy^ Here we explicitly take the domain structure of 
LiHoF4 into account and use a much improved method 
which vastly reduces finite-size corrections. We also use 
a recently introduced cluster algorithms^ In this manner 
we can determine both the critical temperature and crit- 
ical field of the effective Ising model with high precision. 
This precision is high enough that the (considerable) un- 
certainties in the crystal field parameters (described be- 
low) are now the limiting factors controlling uncertain- 
ties in the predicted phase diagram. We will propose a 
simple microwave spectroscopy experiment to eliminate 
these experimental uncertainties. 



II. THE CRYSTAL FIELD HAMILTONIAN 

A single Ho 3+ ion in the crystal LiHoF4 has a par- 
tially filled outermost shell 4/ 10 , and the ground state 
electronic configuration of the Ho 3+ ion is 5 Is as dictated 
by Hund's Rules. The lowest excited electronic configu- 
ration of the ion, lies approximately 7400 K above 
the ground state configuration, as seen in spectroscopic 
experiments on LiHoF4.iS In the range of temperatures 
of interest in this article, any configuration mixing of the 
ground configuration with the excited ones can thus be 
safely neglected. The configuration-mixing due to the 
crystal field is also assumed to be small. 

Considering only the spin-orbit interaction and the 
Hund's Rules, the ground state configuration of the Ho 3+ 
ion in L1H0F4 will be (2J+l=)17-fold degenerate. But 
the interaction of Ho 3+ with the Li + and F~ ions can 
be captured concisely in a crystal-field hamiltonian (Vc) 
that lifts the degeneracy while taking into account the 
symmetry of the crystal. The LiHoF4 crystal has S4 
symmetry, which partially splits the 17-fold degeneracy. 
In S4 symmetry, the states of a configuration with an 
even number of electrons transform according to four one- 
dimensional representations, two of which are related by 
time-reversal symmetry. The ground state of the crystal- 
field hamiltonian is thus a doublet, belonging to the two 
related representations mentioned above, giving rise to a 
non-Kramers degenerate ground state. The crystal field 
hamiltonian depends in a complicated way on the posi- 
tions of the various ions inside a unit cell, but it turns 
out that one can express Vc in terms of the total angu- 
lar momentum (J) operators of the Ho 3+ ions by using 
a set of Steven's equivalent operators and correspond- 
ing phenomenological constants called crystal field pa- 
rameters(CFP) which are determined by fitting to ex- 
perimental spectroscopic and susceptibility data>i2iii(See 



Appendix for details on the crystal field hamiltonian). 



A. Ising system at low temperatures 

Diagonalizing Vc shows that the lowest excited state in 
the spectrum is a singlet, lying ~ 11 K above the ground 
state doublet. At temperatures much lower than this 
gap, only the ground state doublet can be expected to be 
significantly populated, and the low temperature physics 
can be captured by considering a two-state system. The 
two degenerate states (denoted by \a) and \f3}) can be 
chosen such that (a\J z \a) = -{(3\J Z \(3) and (a\J z \[3) = 
0. It turns out that the transverse angular momentum 
operators i x and J v have no non-zero matrix elements in 
the degenerate ground-state subspace. This is the source 
of the strong Ising anisotropy which causes the linear 
susceptibility to vanish in the transverse directions. We 
thus identify the two degenerate states as | f) and | 1), 
and it can be expected that the low temperature physics 
will be described by an effective Ising model with spin-i 
degrees of freedom. 

The S4 symmetry of the crystal defines an easy axis for 
ferromagnetic ordering in the pure LiHoF 4 crystal^ In 
the absence of any externally applied magnetic field, the 
magnetic dipole interaction among the Ho 3+ magnetic 
moments causes them to align along the c-axis of the unit 
cell (the z-direction in this analysis) below 1.53 K and the 
dipolar Ising model serves as an adequate effective model 
for the system. 

The situation becomes more subtle when the crystal is 
subjected to an external magnetic field perpendicular to 
the above mentioned easy axis^ The magnetic moments 
couple to the transverse field through the Zeeman inter- 
action. Restricted within the ground state configuration 
J=8, the Wigner-Eckert theorem yields a Lande g-factor 
()l = \ and the Zeeman term in the hamiltonian can be 
written as 

Hz = -9lVbB ■ J, (2) 

with fiB — 0.6717K/T being the Bohr magneton and 
B = B x e x . Because the Zeeman term has no matrix 
elements within the two dimensional subspace, it cannot 
flip a spin to first order in B. Thus to see any effect 
of the transverse field one must resort to second-order 
perturbation theory. Denoting the singlet excited state 
at A =11 K by I7), one can except to see an effect ~ 

B 2 

{9lHb) 2 -£\{l\J x \oi, ft)\ 2 on the energies of the ground 
states. A naive application of degenerate second-order 
perturbation theory thus suggests that the effect of the 
transverse field should be proportional to B 2 X . 

This perturbation theory scenario breaks down in the 
quantum critical regime as can be seen in a comparison of 
the energies. If we assume that the magnetic moment is 
fully polarized in the transverse direction (( J x ) = J = 8) 
at the quantum critical point of B c x = 4.9T at T=0^ a 
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simple estimate of the Zeeman energy is given by 
E z = 9lHbB c x {J x ) =32.91K, 



(3) 



significantly larger than A. This demonstrates that the 
mixing of the ground-state doublet with all the higher- 
lying states must be considered at large transverse fields 
and second-order perturbation theory is not sufficient to 
incorporate the effect of B in the quantum critical regime. 
We describe below a non-perturbative scheme to capture 
this physics. 



III. MAPPING TO THE ISING SYSTEM 

The magnetic properties of LiHoF4 are determined by 
three kinds of interactions: a long-range magnetic dipole 
interaction among the Ho 3+ magnetic moments, a near- 
neighbor exchange interaction which we assume to be 
small and isotropic, and an isotropic hyperfine interac- 
tion between the electronic and nuclear magnetic mo- 
ments on the same site. Therefore, the complete hamil- 
tonian of a LiHoF4 crystal in a transverse magnetic field 
can be written as 



H = 



i i 



(4) 



C^j contains the position depen- 



where fi, v = x,y,z 

dence of the magnetic dipole interaction so that 



(•5) 



J cx is the dimensionless strength of the antiferromagnetic 
exchange interaction (whose strength is unknown at this 
point), nn signifies that the sum is to be carried out over 
nearest neighbors only and a(=5.175A) is the length of 
the unit cell in the x-y plane. A is the strength of the 
hyperfine interaction(0.039 K) and h is the total angular 
momentum vector of the Ho-nucleus at the i-th site, / = 

7 
2 ' 

To reduce the hamiltonian that acts on a 
(2J+f )x (21+1 )=f 36-dimensional Hilbert space to 
an effective Ising model with spin-i degrees of freedom, 
we neglect the relatively weak hyperfine interaction in a 
first approximation. If A = 0, the only relevant degrees 
of freedom are the electronic angular momenta J, and, 
for J = 8, the single-site Hilbert space is 17-dimensional. 
The transverse field splits the degeneracy in the ground 
state subspace and also mixes the 15 higher-lying crystal- 
field states with the two lowest states. The single-site 
hamiltonian (neglecting the hyperfine interaction) 



is diagonalized numerically for all values of the transverse 
field. For a given B^, let the two lowest states be denoted 
by \a(B x )) and \(3{B X )) and their energies be denoted by 
E a {B x ) and E/3(B X ). Then in a two-dimensional Hilbert 
space which is spanned by \a(B x )) and \f3(B x )} (identi- 
fied as | — >) and | <— ), respectively) Ht can be written 
as 



H T = E CM {B X ) - ^A(B x )cr x , 



(7) 



where E CM {B X ) = \{E a {B x ) + E P (B X )) and A(B X ) = 
Ep(B x ) — E a (B x ). Thus we see that the energy differ- 
ence between the degenerate states caused by the trans- 
verse field can already be interpreted as an effective mag- 
netic field acting on spin-i degrees of freedom at each 
site. Fig. demonstrates how the transverse field con- 
tinuously splits the degeneracy between the two ground 
states. 
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FIG. 1: Lifting of ground state degeneracy by the transverse 
field. 



To generate the complete first-principles hamiltonian 
in Eq. (0J in terms of the Ising spins, one must project 
the ./-operators onto the two dimensional subspace. Rec- 
ognizing the fact that every two dimensional hcrmitian 
matrix can be uniquely expanded in terms of the two di- 
mensional unit matrix and the three Pauli matrices, we 
evaluate the matrix elements of the J-operators within 
the two dimensional subspace for each value of the mag- 
netic field and obtain their representations in terms of 
spin-i operators. Every operator J M (/i=x, y, z) can be 
written as 



J" = C a 



C^{B x )a v . 



(8) 



H T = V c {J)-g L ^BB x r 



(G) 



The Pauli matrices a v requires the choice of a basis | f ) 
and | |) in terms of the states \a) and We choose the 
basis vectors such that | |) = 7^(1^) + exp(z6*) ) and 

| I) = -^(\ a ) — exp(i6*)|/3)). The phase 6 is chosen such 
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that the matrix elements of the operator J z are real. In 
this basis, we have 

Jf = C zz al (9) 

The operators J x and J y are also projected into this ba- 
sis. The following table presents the values of the various 
parameters C M and C^ U1 for B x =4.9 T, which is the ex- 
perimentally observed critical transverse field in the limit 
of vanishing temperature. 








Cx 


3.12 


Cy 


0.01 


Cza: 





Qxx 


0.60 


Cyx 


-0.09 


Czy 





Qxy 


0.01 




1.05 


Czz 


5.14 


Cxz 





Cyz 






TABLE I: The strengths of various coefficients C MI / for B^=4.9 
T. 



of ~ 1%, the effective Ising model for L1H0F4 can now 
be written as 

iflsing = ^ y E C M,i(B x ) — -^A(B X ) ^ of 

i i 

+ \{g^ B c zz f^Y.«n- (11) 

i,nn 

The Ising hamiltonian of Eq. (|llf> acts on effective spin- 
\ objects("pseudospins") located at the lattice sites of 
Ho 3+ ions in the pure crystal LiHoF4, and the interac- 
tion strengths and the effective transverse field depend 
on the physical transverse field B x implicitly through the 
parameters C zz and A. The parameter Eqm provides a 
B-j-dependent zero of energy and will be ignored for the 
purpose of computing the magnetic phase diagram. 



We replace the J-operators by their corresponding rep- 
resentations in the two-dimensional basis in the hamilto- 
nian in Eq. (@J (The hyperfine interaction is neglected 
at this stage). This leads to an extremely complicated 
hamiltonian that acts on the Ising subspace. The projec- 
tion generates various new kinds of interactions among 
the effective Ising spins. For example, the term 

Hzx = -^{9l^b) 2 ^ CjjJfJj 

in Eq. |0J gives rise to a set of complicated interactions 
in the Ising subspace that are written as 

H zx = \(9LHB? Y, £™{C zz ot)(C x + C xx a* + C xy af) 

i jjiO 

+\{gLHB) 2 Y. c ™ CzzC *y^l ( 10 ) 

In Eq. JU3J, the quantities C zz etc. depend on the trans- 
verse magnetic field, B x . The strength of the effective 
interactions generated by the Ising projection is deter- 
mined not only by the parameters C M „, but also by the 
parameters C"? , that depend on the inter-spin vector fij . 
Thus we can compare the strengths of the various effec- 
tive interactions, and we find that the largest effective 
interaction is J*Jj oc (C zz ) 2 a z a z . This interaction is 
larger than the other interaction terms by two orders of 
magnitude for the entire range of magnetic fields in ques- 
tion (except for some constant terms that just serves to 
define a new zero of energy). Therefore, to an accuracy 



IV. COMPUTING THE PHASE DIAGRAM 

Eq. (|llf) is a particular example of a general class of 
hamiltonians in which the various terms in the hamilto- 
nian do not commute with each other, and in the interest- 
ing parameter regime around the quantum critical point, 
they are comparable in strength. Therefore, quantum 
fluctuations can become strong enough in the system to 
destroy long-range order even at T=0. An Ising model 
of spin-i objects on a three dimensional lattice placed 
in a magnetic field transverse to the Ising direction is 
a very good prototype of such a hamiltonian. We uti- 
lize Eq. (|llf> as a starting point to compute the phase 
diagram of LiHoF4 in various ways, notably using mean- 
field theory and quantum Monte Carlo simulations. As a 
spin-off, we are able to calculate the strength of the ex- 
change interaction J cx . Finally, we compare our results 
with existing experimental dataj2> which were obtained 
from magnetic susceptibility measurements. 

A. The effective hamiltonian 

An unusual feature of the Ising hamiltonian in Eq. I|ll|) 
is that the strength of the dipole interaction itself seems 
to depend explicitly on the strength of the physical trans- 
verse field through the parameter G zx . Ordinarily, the 
interaction term in an Ising model is free of any depen- 
dence on external magnetic fields. The popular hamil- 
tonian governing the quantum Ising model is generically 
written as 

where J^- is the pairwise interaction term (dipole interac- 
tion, near- neighbor exchange ..) which usually depends 
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on the difference in spatial positions of the spins, and h x 
is the magnetic field that acts as the source of the quan- 
tum fluctuations. With some appropriate definitions, we 
are able to define an effective model H g that is of the 
same general form as Eq. I|12|l . 

The reduction to a two dimensional Hilbcrt space from 
a 17-dimensional Hilbert space gives rise to renormalized 
Lande g-factor g||, defined as 

.9,1 = 2g L (a(B x = 0\J z )\a(B x = 0)> = 13.8. (13) 

The factor 2g^C zz {B x ) can be interpreted as a renormal- 
ization of the individual magnetic moments due to the 
presence of the transverse magnetic field. To extract this 
renormalization factor, we define a dimensionless ratio 
e(B x ) such that 



e(B x ) 



2g L C zz (B x ) 
3|| 



(14) 



In the classical limit, B x — 0, e(B x ) = 1 by definition. 
This makes it possible to define an effective hamiltonian 
H e ff, which is essentially identical to Eq. (|12[) . 



[e(B x )] 2 H eS , 



(15) 



where 



Htff — — 



2^ 



9\\^b 



+ 



1 

2 ^ 



JfiZ 



A(B X ) 



2 



(16) 



We define every length in the system in units of a = 



5.175^4, such that 



— . Then the spatial depen- 



dence of the dipole interactions can be expressed in 
the terms of dimensionless quantities Ly 1 ', defined as 

ij a- 



-. Replacing the constants in Eq. I|16[l . we 

see that (^) 2 4, - 0.214 K, and (^) = 4.635 
K/T. These constants specific to LiHoF4 are substituted 
in Eq. (|16f) . and we obtain 



H, 



cir 



zz, 
'ij ' 



i x 0.214 J**Pi<>j 

i,nn 

4.635 J B£ f 5>?. 



(17) 
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FIG. 2: The effective magnetic field. 



We note that the dependence of B* s on B x in Fig. (J2J 
is indeed quadratic for small values of B x , as dictated 
by second-order perturbation theory but the depen- 
dence becomes increasingly weaker as B x increases. This 
demonstrates that in the quantum critical regime of large 
transverse fields and vanishing temperature, perturba- 
tion theory arguments are not adequate. 

At this point, a subtlety in the definition of "temper- 
ature" needs to be mentioned. All our further calcula- 
tions use H e ff as in Eq. i|17|) . but the physical hamilto- 
nian is -ffismg, which modifies the energy scale by a factor 
[e(B x )) 2 (See Eq. ifPo|l % ). In both mean-field theory and 
quantum Monte Carlo calculations, the temperature en- 
ters only in the definition of the Boltzmann weight of a 
configuration or a state- vector in the form (ke = l) 



p = exp 
= exp 
= exp 



His 



T 

[e{B x )fH eS 



T 



H, 



eff 



T. 



c rr 



(19) 



where the physical temperature T is related to the effec- 
tive temperature T e g by 

T = [e(B x )] 2 T cS . (20) 



B. A possible three-state description of the 
effective Ising model 



A comparison of Eq. (|16fl and Eq. I|17l) immediately yields 
the crucial correspondence between the effective trans- 
verse magnetic field B* s (expressed in Tesla) and the 
physical transverse field B x (also expressed in Tesla). 



2 x 4.635 x [e(B x )} 2 



(18) 



In this section, we digress from the discussion a bit to 
show that it is also possible to capture the degeneracy 
splitting A(B X ) using a simple 3-state model of the sys- 
tem, but it leaves out the effect of the matrix elements 
of the ground state doublet with the higher-lying states. 
If we consider the states |a),|/3) and I7) at B x = as 
the basis states, then the crystal field hamiltonian Vc is 
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simply diagonal in this basis, and looks like 



V C 



( 0.0 
0.0 
v OA 



(21) 



where A = 10.8K. 

The operator J x has only two non-zero matrix ele- 
ments restricted within the three states in the absence 
of a transverse field. By choosing the phase of the wave- 
functions properly, it is possible to make both of them 
real. These matrix elements are 



(a|J*| 7 > = </3|J*| 7 )=p, 



(22) 



where p — 2.4. Then the Hamiltonian in Eq. 10 can be 
written in the three-state basis as 



1 -gLHBB x p 

-g L p B B xP | . (2:-!) 

\ -9lPbB x p -QlUbBxP A 



This hamiltonian can be exactly diagonalized, and defin- 
ing B x = ghPBB x p, the eigenvalues are 



E 



1 



8B 2 



2 2 

E 1 = 

E 2 = ^ + ^A 2 +8^. 

The eigenstates can be schematically written as 
1 



(24) 



|V>o) 



V2 
1 

72 



(|a> + |/3»+e |7> 
(l«) ^ 1/3)) 



(25) 



1^2) = \j)+e 2 (\a) + \P)), 

where the symbols Eq and e 2 signify small admixtures. 

It is readily seen that the transverse field splits the de- 
generacy by an amount A(B X ) = E\ — Eq. The gap is 
quadratic in B x when B x <C A, but it is, in fact, linear 
when B x ^> A. However, we see that the lowest excited 
state has an eigenvalue of zero for all B x , unlike the real 
system, where it bends downwards with increasing B x . 
Thus the three state model cannot capture the repul- 
sion from the higher-lying states. Even then, it explicitly 
shows that the physics near the quantum critical point is 
ill-described by second-order perturbation theory. 



C. Domains in LiHoF4 and the mean-field solution 

The magnetic dipole interaction is long-range, falling 
off as \, where r is the distance between two spins. A 
dipole-coupled Ising ferromagnet is expected to be well 
explained by mean- field theory for d > 3. As we shall see, 



mean-field theory does capture the qualitative behavior 
in the phase diagram, but is not very reliable for making 
quantitative predictions. 

Our starting point is Eq. (|17f) . The individual spin-i 
degrees of freedom sit on a tetragonal unit cell of dimen- 
sions (a, a, c) = (1, 1, 2.077) in units of a = 5.175A Each 
unit cell has four spins, at the positions we denote as 
[(0,0, 0), (0, §, f),(§, § , -f ) and (f ,0, -f )]. All the po- 
sitions in the unit cell are measured from the Bravais 
lattice vector denoted by R + (0, 0, 0). 

The long-range nature and angular dependence of the 
dipole interaction create a complex model in which the 
actual shape and lattice of the sample influences the 
ground state. Spins that lie in the same ab-plane want to 
antialign, while spins along the same c-axis want to or- 
der ferromagnetically. This angular dependence tends to 
favor ferromagnetism in long thin samples, while it pro- 
hibits a ferromagnetic state in spherical samples. Lut- 
tinger and Tisza— found, for example, a ferromagnetic 
ground state for the face-centered lattice in the shape 
of a prolate spheroidal sample of axis ratio larger than 
6. However, in experiments on LiHoF4 spherical, 2 cubic 4 
and rectangular 1 : shapes have been used with no apparent 
dependence of the results on the shape of the sample. 

The reason for this sample-shape independence lies in 
the domain structure of LiHoF4. Experimentally, there 
is evidenc o 13 i 14 that LiHoF4 forms long needle-shaped 
domains, with the long direction being along the c-axis. 
Since this is an Ising system, the magnetization alter- 
nates in sign between adjacent domains, which minimizes 
the global energy stored in the macroscopic fields outside 
the sample. Thus it would be incorrect, in the context 
of mean-field theory, to assume that the magnetization 
is uniform over the entire sample, since in fact, it is only 
uniform over a single domain. 

To incorporate the effect of domains, we assume that 
an imaginary sphere sits deep inside a single domain^ A 
single spin at the center of this sphere now experiences an 
effective field from the orientation of the magnetic dipoles 
both inside and outside the sphere. In the ferromagnet- 
ically ordered state, the magnetization of the domain in 
question can be assumed to be uniform, denoted now by 
M 2 . The magnetic field acting on the domain as a whole 
is the external field i? 2 xt , and the susceptibility \ °f the 
domain can be found from the relation 



M z = X Bl 



(26) 



On the other hand, if we consider the small imaginary 
sphere deep inside the domain, the susceptibility Xsph of 
the sphere can be found from the equation relating the 
magnetization inside the sphere (since the sphere is a 
part of the domain, the magnetization inside it is M z ) 
and the magnetic field acting on the sphere, B z ph , 



M 2 



Xsph-S S ph 



(27) 



We treat the magnetic field produced by the spins out- 
side the sphere using mean-field theory. The spins at the 
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surface of the hollow sphere produces a polarization field 
proportional to M z , while the spins at the outer surface 
of the domain produces a demagnetization field, also pro- 
portional to M z , that is characterized by a demagneti- 
zation factor N z , dependent on the shape of the domain. 
Therefore, we can write down the general relation^ 



mean-field given by 



B. 



ph 



B z + ^-M z - N Z M Z 



N z 



4tt 



(28) 



and sub- 



For a thin needle-shaped domain, 
stituting Eq. in Eqs. JSjjJ and leads to the fol- 
lowing relation 



1 



X = 



Xsph 



4tt 
3 



(29) 



Eq. (|29|l above is the magnetic version of the Clausius- 
Mosotti relation, 16 which relates the macroscopic electric 
susceptibility of a system with the microscopic molecular 
polarizability. In fact, it is in some ways easier to derive 
from the analogy with electric dipoles. In that case, the 
field produced by the charge density at the surface of the 
spherical cavity is ^-M z while the small random patches 
of charge where the domains reach the surface of the 
sample produce negligible field. 

We emphasize that this relation is obtained from a 
single assumption that the spins outside the sphere are 
treated in mean-field fashion. The microscopic variable, 
Xsph, can be calculated approximately, using mean-field 
theory, or exactly, using quantum Monte Carlo(QMC) 
simulations over a sphere. The next part of this section 
explores how the domain structure can be incorporated 
in the mean-field scenario, while the following section 
contains the results from QMC simulations. 

At this point, a comment on the effect of the trans- 
verse magnetic field on the domain formation is neces- 
sary. From the mapping to the Ising model, it is seen that 
the contribution of the transverse dipole interaction, Lf^ 
is negligible compared to the longitudinal dipole inter- 
action, Lfj. Therefore, we may assume that the trans- 
verse magnetic field polarizes the Ising spins along the 
x-dircction uniformly throughout the sample. In other 
words, the x-component of the magnetization, M x , is 
unaware of the domain structure formation in LiHoF4. 



D. Mean-Field theory 

Let us now calculate the mean-field critical temper- 
ature for a long, needle-shaped domain. We therefore 
place the sphere deep inside a single domain. The mean- 
field critical temperature is given by the local field, which 
is given by the sum over all dipoles within the sphere. In 
the presence of domains, this field must be augmented 
by the field acting on the sphere, B z ph . 

With no domains a single spin at the center of the 
sphere therefore experiences an effective longitudinal 



B 



MF 




(30) 



where m z is the magnetic moment of a single dipole, 

m-=(»)<o«>. (3D 

The value of the lattice sum in Eq. I|30|l can be easily 
obtained by summing over larger and larger spheres, and 
the convergence is found to be very rapid. The value for 
LiHoF4 is 



sphere 

E 



= 3.205. 



(32) 



In order to take the domain into account we combine 
Eqs. l(2"5)l and Eq. l|3*U|l . Since the macroscopic magneti- 
zation is given by 



M z = N n m 2 



(33) 



where Nq is the number density of dipoles, we simply 
have to augment the value of the lattice sum in Eq. I|32|) 
by the value 4fiV a 3 . Since N = V 4 V and c/a = 
2.077, it follows that ^-N a 3 = 8.067. 

This leads to the value of the critical temperature at 
the classical limit. 



T c = 0.214K x (3.205 + 8.067) w 2.41K. 



(34) 



On the other hand, the new critical effective field at 
the quantum limit is given by 



BX 2.41K 
cff ' c 4.635K/T 



0.52T, 



(35) 



which corresponds to a physical critical transverse field 
ofB x , c «4.11T. 

We note that the critical temperature at the classical 
limit, 2.41 K, is 57% larger than the experimentally ob- 
served value of 1.53 K. This behavior is encouraging, in 
the sense that mean-field theory should always overesti- 
mate the ordered region. As far as the quantum limit is 
concerned, it may seem that the mean-field estimate is 
wrong, since 4.11 T is about 16% smaller than the ex- 
perimentally observed value of 4.9 T. However, we have 
entirely neglected the on-site hyperfine interaction be- 
tween the electronic cloud and the Ho-nucleus. At low 
temperatures, this interaction affects the phase bound- 
ary significantly- Therefore, at this point, we cannot yet 
compare our results for the quantum limit with experi- 
ment. 

As a consistency check to our mean-field treatment of 
the domain structure, we perform the lattice sum directly 
over a long cylinder. The summation is done by consid- 
ering a series of long cylinders with increasing base area. 
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For each cylinder with a fixed base area, the length is 
increased until convergence is obtained. The series con- 
verges quite rapidly and as a final result we have 

cylinder / ~ 2 q z2 \ 

J2 \%-JZ- \ = H.272 



j¥0 



1 j 



3.205 + 8.067, (36) 



in perfect agreement with the mean-field approach de- 
scribed above. We note that the mean-field theory cal- 
culation in Ref . [2] does not take into account the physics 
of domain formation but rather simply re-scales the ef- 
fective couplings to force a fit to the observed Tc and 
B x n- 



V. QUANTUM MONTE CARLO SIMULATIONS 
ON LIHOF 4 

We perform extensive QMC simulations of the dipole- 
coupled quantum Ising model, described by the hamilto- 
nian in Eq. (|17l) . We use a recently introduced^ stochas- 
tic series expansion (SSE) cluster quantum Monte Carlo 
method for which computational time scales as Nln(N), 
where the number of spins is given by N. In traditional 
single-spin- flip simulations of long-range models the com- 
putational time typically scales as N 2 , due to the sum- 
mation over interactions between all pairs of spins. The 
improved scaling enables us to increase the precision and 
reach large system sizes. 

Our starting point is again Eq. I|29|l . which relates a 
microscopic quantity, Xsph, to a macroscopic quantity x- 
From Eq. I|29|l . we find that x can diverge, even when 
-Bf xt vanishes, when 



Xsph 



_3_ 

47r' 



(37) 



We proceed to evaluate Xsph exactly using QMC simula- 
tions as a function of T e g and -Bjfff- At the point where 
the condition in Eq. (|37fl is met, x, the macroscopic sus- 
ceptibility, diverges and the system is critical. 

Truncating the sum in the hamiltonian at the bound- 
ary of the sphere and treating the dipoles outside the 
sphere in mean-field fashion can be considered a bound- 
ary condition in the Monte Carlo simulation. This so- 
called reaction-field methodic is not the only way to in- 
corporate the long-range interaction in an efficient way. 
An alternative is to consider a large number of periodic 
images of the simulation volume. Performing the nec- 
essary sums over all dipoles in the simulation volume 
and the image volumes is very time consuming, but the 
sums can often be efficiently performed using the Ewald 
summation technique^ In this work we have not at- 
tempted to compare the relative advantages of the meth- 
ods, but we have found the reaction-field method very 
convenient. An earlier comparison^ found the Ewald 
summation method also quite reliable. 



Operationally, the susceptibility Xsph is defined as the 
spin-spin correlation function in imaginary time, aver- 
aged over all spins inside the sphere. 

X^ = ^Y,[dr{at{r)a^)). (38) 



The prefactor a is given by 



(39) 



The imaginary time integral in Eq. (|38|l can be evaluated 
directly by the SSE method^ The condition for critical- 
ity can be now be written as 



4?riVo \g\\HB 



= 0.579 KT 



(40) 



The simulation is done over a sequence of spheres with 
increasing radius. We find that the critical curve given 
by the above condition converges fairly quickly as the 
radius is increased. In the QMC calculation the fluctua- 
tions within the sphere are fully included, and we stress 
that although we use a mean-field result for the part of 
the domain exterior to the sphere, the final result has 
converged and should not, in any way, be considered a 
mean- field result. One could also argue that the field 
originating from dipoles in other domains may affect the 
critical temperature in LiHoF4. However, the material 
forms domains in order to minimize the energy density 
of the magnetic field, so this field should be very small. 
The same arguments hold also for the case of a transverse 
field, and the same condition is used in the presence of 
a transverse field. Fig. © shows the phase diagram for 
the mean-field solution as well as for the Monte Carlo 
simulation. 

The results of Fig. © seem reasonable in the sense 
that the QMC solution yields a lower critical temperature 
than the mean-field solution. The mean-field solution is 
within 20% of the QMC solution, and it does therefore 
not suffice to describe the system at a more precise quan- 
titative level, even for the present case of a long-range 
interaction in three dimensions. The finite-size effects 
in Fig. are quite small already for moderate system 
sizes and the result for the largest system size (N=295 
dipoles) is T c = 2.03 K at 13^=0 and {B* s ) c = 0.47 T 
at the quantum limit. The critical temperature is still 
significantly above the experimental value of 1.53 K, but 
this discrepancy may be attributed to an additional ex- 
change interaction, which we discuss shortly. The effec- 
tive transverse field of 0.47 T corresponds to a physical 
transverse field of approximately 3.77 T. But, as stated 
earlier, comparison of the critical field with experimen- 
tal value at the quantum limit is possible only when the 
hyperfine interaction is included in the model. 
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FIG. 3: The phase diagram as a function of effective tem- 
perature and effective magnetic field obtained from quantum 
Monte Carlo simulations on different lattice sizes. Upper 
curve is the mean field solution. 



A. Exchange interaction 

A major question at this point is why the critical tem- 
perature of the model is about half a Kelvin higher than 
the experimental value. One factor that lowers the crit- 
ical temperature is the antiferromagnetic Heisenberg ex- 
change interaction. Since there has been no direct ob- 
servation of the strength of the interaction, we treat 
it as a free parameter, J ex , as in Eq. (|17|l . First we 
consider the mean-field solution with the domain struc- 
ture. If we assume an exchange interaction of strength 
J cx = 1.03 the mean- field critical temperature is obtained 
as T c = 0.2140K x (11.271 - 4 x 1.03) = 1.53 K, since 
every spin has four nearest neighbors. 

Next we perform Monte Carlo simulations, where we 
tune the value of the exchange interaction to obtain a 
critical temperature of 1.53 K. This is done for a se- 
quence of spheres with increasing radius, and the required 
exchange parameters are shown in Fig. The corre- 
sponding phase diagram is shown in Fig. JSJ). There is 
still a strong (and non-monotone) size dependence in the 
exchange parameter due to fluctuations in the number of 
surface bonds as the radius of the sphere increases, but 
the phase diagram shows very little size dependence. The 
exchange parameter itself has not converged for the sys- 
tem sizes we study, but the results for the largest system 
sizes averages to about 0.75. This is about 25% smaller 
than the mean-field value, which seems reasonable. A 
value of J cx = 0.75 corresponds to an antiferromagnetic 
exchange interaction of strength 0.16 K between neigh- 
boring spins. This can be compared to the dipolar inter- 
action between nearest neighbors, which is ferromagnetic 
and of strength 0.31 K, or about twice the exchange in- 
teraction. If we instead sum over all bonds connected to 
a given site in an ordered domain, the dipolar interac- 
tion is of strength 0.214 x 11.271K = 2.4 K, while the 
exchange energy is of order 0.214 x 4 x 0.75K = 0.64 



K. It is difficult to know if this is a reasonable value for 
the exchange, and a more stringent test comes when the 
whole phase diagram in can be compared to the experi- 
mental result. In the next section, we consider including 
the hyperfine interaction, which enables us to compare 
the phase diagram to experimental data. 

The reason for the strong finite-size effect in the ex- 
change parameter is that the exchange interaction for all 
broken bonds at the surface of the sphere is neglected in 
this calculation. If a spin is located close to the boundary 
of the sphere then one, two or three of its four nearest 
neighbors may be located outside the sphere. Even for 
the largest system size (N=3491) in Fig. (@|, only 80% 
of the spins have all four nearest neighbors inside the 
sphere. The fraction of spins that have only one, two or 
three nearest neighbors inside the sphere fluctuates very 
rapidly as the size of the sphere is increased. However, 
this is a boundary effect and should disappear as the sys- 
tem size is increased further. 
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FIG. 4: Exchange energy needed to tune the critical temper- 
ature. 
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FIG. 5: Critical temperature as a function of effective trans- 
verse field, with an exchange interaction 
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VI. THE HYPERFINE INTERACTION 

The dynamics of LiHoF4 at low temperatures is com- 
plicated by the hyperfine interaction between the elec- 
trons in the Ho 3+ ion and the Ho-nucleus^ This inter- 
action is assumed to be on-site, and its strength is char- 
acterized by the hyperfine constant, A, which is equal to 
0.039 K for LiHoF 4 i 21 i 22 i 23 At low temperatures, the hy- 
perfine interaction prefers ordering and the final phase 
boundary is a result of the competition between the 
transverse field , which tries to destroy the ferromagnetic 
ordering, and the hyperfine interaction. 

Consider the truncated hamiltonian in Eq. ©. Since 
the value of the total nuclear angular momentum I in 
LiHoF4 is |, the actual Hilbert space of the electron- 
nucleus system for a given site is, in reality, 17 x 8 = 136- 
dimensional. With the inclusion of the hyperfine cou- 
pling, the truncated hamiltonian in Eq. © should be 
written as 

ffiw = (Vc(J) - 9lVbB x J x ^ ®l N + A J2 ® -P*, 

fi— x,y,z 

. ( 41 ) 

where Iat is the eight dimensional unit matrix in the 
nuclear sector of the Hilbert space. 

The hamiltonian in Eq. I|41[) is a 136-dimensional ma- 
trix, and in the absence of the hyperfine interaction, each 
electronic crystal-field wavefunction is 8-fold degenerate. 
The hyperfine interaction breaks this degeneracy to the 
order of the strength of A. We adopt the viewpoint that 
the transverse magnetic field is renormalized because of 
the hyperfine interaction, in a fashion controlled by the 
temperature. This renormalization can be extracted by 
matching the longitudinal susceptibility of a single Ho 3+ 
ion without the presence of the hyperfine interaction as 
described below. 

The susceptibility, \zz , for a single ion is easy to define 
in the 136-dimensional Hilbert space. We diagonalizc 
the Hamiltonians in Eq. (gTJ with A=0 and A=0.039 K. 
Let the eigenstates in the 136-dimensional Hilbert space 
be denoted by [tp^ > when A=0 and [ip^ > when A is 
nonzero, and the corresponding energy eigenvalues as E^ 
and E^ n , m=l . . . 136 in both cases. For a temperature 
T = 4, the longitudinal susceptibility can be defined as 



z 



1 m,n— 1...136 

-2 A \{r m \J z ®lNm\ 2 e-^ 



E 



m,n—l. . . 136 



EL - El 



(42) 



where i=0,l. The primed sum is over states degener- 
ate in energy, while the double-primed sum is over states 
non-degenerate in energy. Fig. JfjJ shows how the sus- 
ceptibilities at the typical temperature T=0.444 K differ 
for a range of values of the transverse field depending on 
whether the hyperfine interaction is turned on or not. 




FIG. 6: Comparison of the single-ion longitudinal suscepti- 
bilies with and without the hyperfine interaction. We make 
the approximation of assuming that the increase in longitu- 
dinal susceptibility caused by the nuclei can be modelled by 
simply renormalizing the transverse field downwards by the 
appropriate amount. 



From the plot of the susceptibilities, it is clear that 
the hyperfine interaction renormalizcs the magnetic field. 
This effect is temperature-dependent, and the shift in 
susceptibility decreases as the temperature is increased. 
For instance, the Ising mapping and quantum Monte 
Carlo simulations yield a typical point on the phase 
boundary as (T C ,B X)C ) = (0.444K, 2.949T). In Fig. © 
we see that xf (A = 0, T = 0.444K, B x = 2.949T) equals 
Xf (A = 0.039K,T = 0.444K,B X = 3.282T). Thus, we 
conclude that the critical transverse field of 2.949 T is 
renormalized to the critical transverse field of 3.282 T 
in the presence of the hyperfine interaction. This sim- 
ple renormalization program can be carried out for all 
the points on the quantum Monte Carlo phase bound- 
ary. This procedure yields a phase diagram that shows 
how the hyperfine interaction affects the phase boundary 
significantly in the quantum regime. 



A. Magnetic phase diagram from quantum Monte 
Carlo and hyperfine interaction 

Fig. shows a comparison of the theoretically ob- 
tained phase diagram with the experimental phase dia- 
gram obtained through susceptibility experiments. 2 We 
get a quantum critical point of B x .c — 4.66T, which is 
approximately 6% smaller than the experimental value. 
However, serious deviations also occur at temperatures 
higher than T ~ IK, where the effect of the nuclear in- 
teraction is negligible. This leads us to conclude that 
the deviation really stems from the mapping to the Ising 
model. The Ising model mapping is strongly determined 
by the strength of the degeneracy splitting in the ground- 
state doublet as the transverse magnetic field is turned 
on, which in turn strongly depends upon the values of the 
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crystal field parameters. In a subsequent section, the un- 
certainties in the crystal field parameters and their effect 
on the physics are discussed. 



A± has any non-zero matrix element within the same de- 
generate electronic multiplet. Thus An drops out of the 
physics in first-order perturbation theory. 



B. Hyperfine interaction and the Ising model 

The actual computation of the longitudinal susceptibil- 
ity of the ion-nucleus system in the presence of a hyper- 
fine interaction requires the knowledge of all the eigen- 
states of the 136-dimensional Hilbert space. However, the 
physics of why the hyperfine interaction enhances the lon- 
gitudinal susceptibility can be understood by considering 
a toy system consisting of a single spin-i electron coupled 
to a spin-i nucleus through a hyperfine interaction. It 
shall be seen that the enhancement of susceptibility is re- 
ally a subtle effect captured in second-order perturbation 
theory. 

Let us consider a single spin-i electron placed in a 
transverse magnetic field h x (we assume that the nucleus 
does not couple to the magnetic field: this is a reasonable 
assumption since the nuclear g-factor in LiHoF4 is ap- 
proximately 1000 times smaller than the electronic Lande 
g-factor). Then the hamiltonian is 



H 



1 



N- 



(43) 



The subscripts e and N refer to the electronic and nu- 
clear degrees of freedom, respectively. The Hilbert space 
is ((2x1 + 1)(2 x i + 1) =)4-dimensional, and the 
transverse field splits the Hilbert space in two multi- 
plets, each of which is twofold degenerate. The ground 
electronic multiplet consists of the states | — *) e <8 | 1)n 
and | — *) e (g) | 1)n, while the excited electronic multiplet 
consists of the states | «— ) e <8> | 1}n and | «— ) e <8 | |)jv, 
the difference in energy between the two multiplets be- 
ing 2h x . The longitudinal susceptibility can be easily 
evaluated using the expression in Eq. (|42|l (with suitable 
modifications for the toy system), and we find 



Xzz — , 

h :r 



(44) 



at T = 0. 

If the hyperfine interaction is turned on, the degenera- 
cies between the states in the nuclear sector is lifted in 
each multiplet, and this changes the susceptibilities. Let 
the perturbing hyperfine interaction be written as 



Vhyp — All 



>N 



A 



±<J e <J N - 



(45) 



Even though the hyperfine strengths are isotropic in the 
physical system, introducing anisotropy helps us to un- 
derstand the roles played by the longitudinal and trans- 
verse components of the hyperfine interaction transpar- 
ently. The special case of isotropy, An = A±, will be 
considered at the end. 

In first order degenerate perturbation theory, the elec- 
tronic states are all polarized in the x-direction, and only 
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FIG. 7: The complete phase diagram of LIH0F4. Experimen- 
tal data is from Ref.0. 
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FIG. 8: Ion energy spectrum (for a spin 1/2 nucleus) with 
the hyperfine interaction treated in first order perturbation 
theory 

In Fig. JSJ), the states in the four dimensional Hilbert 

space are = | ->) e ® | <-)jy, IV'a) = I "Oe ® I ~*)n, 
1^3) = I <-}<= ® I -+)n and \il>%) = I Oe <8 I HiV- The 
zero temperature susceptibility is now given by 



1 



Xz 



h x + I Aj_ 



(46) 



Thus it is obvious that the longitudinal susceptibility 
always decreases in first-order perturbation theory. The 
transverse component of the hyperfine interaction always 
acts as an extra transverse field, which implies that the 
hyperfine interaction should lower the critical field, at 
least in first-order perturbation theory 

However, if we consider second-order perturbation the- 
ory, we can show that, when the hyperfine interaction is 
isotropic, the net effect is actually an increase in the lon- 
gitudinal susceptibility at T = 0. In second-order, only 
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An contributes to further splitting of the energy spec- 
trum. More importantly, the ground-state | — >) e ® | «— )jv 
is mixed with the state | «— ) e ® | — >)at, and the state 

I — *)e ®\ —>-)n is mixed with the state | <— ) e ®\<—)n 

An 

both with an amplitude given by £ = t-"-. Thus, up to 
second-order perturbation theory we have the following 
spectrum: 



1^4) 



An 



An 



1^4) 



1^3°) + 



(47) 



The energies of the states up to second-order in pertur- 
bation theory are given by 



£1 



-ft, - A ± - 



E 2 = -h x + A ± - 



E, 



-ft, - A i 



E 4 = +h x +A± 



2h x 

3. 

2h x 

if. 

2^ 

iL 

2ft. ' 



(48) 



This gives rise to a non-zero matrix element of erf be- 
tween \ipi) and \1jj2), and this term serves to cancel the 
decrease in susceptibility because of first-order contribu- 
tion from A±. The exact expression for the longitudinal 
susceptibility at T — is given by 



2 \2 



q-(^) 2 ) 

2h x 



At 



A 1 



2h x A± 



(49) 



The denominator in the first term in the RHS of Eq. H49|) 
can be easily expanded in terms of the ratios -r^ and 

A 2 

2$?, and the susceptibility expression turns out to be 
(C > A\\,A_l) 



1 



A\ 



+ — - 

h x A± 



A 
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'ftl 



(50) 



The crux of the increase in susceptibility lies in the 
fact that the term proportional to is positive when 

A\\ > A±. The increase in susceptibility is not a generic 
feature of the ion-nuclear interaction for all values of Au 
and Aj_, but is present when the interaction is isotropic, 
as in LiHoF4. The transverse component always drives 



the susceptibility lower, but the longitudinal component 
can compete and win for a certain range of values, and the 
change in susceptibility is quadratic in Ay . We emphasize 
that this positive contribution arises from the non-zero 
matrix element between two states that were degenerate 
in the absence of the hyperfine interaction. 

This schematic problem of a spin-i electron coupled to 
a spin-i nucleus through a hyperfine interaction serves 
nicely to illustrate the physics behind an increase in the 
longitudinal susceptibility of the electrons due to the hy- 
perfine interaction. In LiHoF4, every electronic state 
is split into eight nuclear states when the hyperfine in- 
teraction is turned on. There is no matrix element of 
J z within the same electronic multiplet when Am = 0. 
When An 7^ 0, a non-zero matrix element of J z devel- 
ops between two lowest nuclear states in the ground elec- 
tronic multiplet, the square of the matrix element varying 
quadratically as Ay . Exactly as outlined in the schematic 
problem, this matrix element causes an increase in the 
longitudinal susceptibility of the electron. 



C. Projecting the hyperfine interaction 

The hyperfine interaction involves the electronic angu- 
lar momentum J M . Following the mapping to the Ising 
subspace in the previous sections, we can map these de- 
grees of freedom to effective Ising degrees of freedom. 
This gives rise to a hyperfine interaction term that op- 
erates in an effectively 16(=2x 8)-dimensional Hilbcrt 
space. However, this gives rise to a hyperfine interaction 
whose strength depends on the magnetic field. Using the 
mapping in Eq. Q, Eq. JSJ) and the hyperfine interac- 
tion in Eq. I|41|l . the effective single-site hamiltonian now 
becomes 



H 



a x e ®lN + AC af} {B x )o- a ®lf. (51) 



Using the projected hamiltonian in Eq. I)5ip. we can 
naturally understand why the longitudinal component of 
the hyperfine interaction wins over the transverse one. 
The longitudinal component is proportional to C ZZ (B X ), 
which is significantly larger than all the other compo- 
nents at the magnetic field regime of interest. The 
analysis of the previous section remains valid in a 16- 
dimensional Hilbert space, and the susceptibility is al- 
ways enhanced due to the contribution from the lon- 
gitudinal component of the hyperfine interaction A11 = 
AC ZZ (B X ). 



D. Uncertainties in the Crystal Field Parameters 

The quantitative details of the LiHoF 4 depend sensi- 
tively on the values of the various crystal field parame- 
ters (CFP) used in constructing the crystal field hamil- 
tonian (Please see Appendix for details). For example, 
we have found that the value of the critical transverse 
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field, B Xy c, varies by as large as 25% when the CFP 
is varied by 10%. However, the parameters cannot 
be determined directly by experiments, but are used as 
fitting parameters to fit theoretical calculations to ex- 
perimental data. There have been attempts to deter- 
mine the CFP-s by fitting to spectroscopic datai^ or to 
susceptibility measurements. 11 However, we have found 
that the results of theoretical calculations become more 
and more sensitive to the values of CFP-s as the trans- 
verse field is increased. As we have shown in earlier sec- 
tions, the most important physical quantity that deter- 
mines the phase diagram is the splitting between the two 
lowest states, A(B X ), that smoothly and monotonically 
increases with the transverse field. Thus we propose a 
spectroscopic experiment in the presence of a transverse 
field to determine A(B X ). However, the quantum regime 
of the phase diagram is complicated due to the presence 
of hyperfine interaction. The effect of the hyperfine in- 
teraction is also extremely sensitive to the values of the 
CFP-s, and this effect is very difficult to calculate theo- 
retically. Thus we propose the experiment be carried out 
in a regime of the phase diagram where the quantum fluc- 
tuations due to the hyperfine interaction are negligible, 
but the fluctuations caused by the transverse magnetic 
field are still significant. Thus we propose an experi- 
ment to determine A(B X ) in the magnetic field range 2.0 
- 3.0 T. In this regime, our theoretical calculations show 
that the splitting energy varies between ~ 1.6 K - 3.0 
K. This corresponds to a microwave frequency range of 
^20-70 GHz. A spectroscopic experiment that deter- 
mines A{B X ) accurately is the most important ingredient 
needed to determine the phase diagram accurately. Even 
though a single spectroscopic experiment is not enough 
to determine the values of all the crystal field parameters 
uniquely, we have found that the only other important 
parameter in the effective Ising hamiltonian, C ZZ (B X ), is 
extremely robust even to large changes in the crystal field 
parameters. Thus, an accurate determination of A(B X ) 
is enough to compute the phase diagram from the effec- 
tive Ising model. 



VII. SUMMARY AND CONCLUSIONS 

It has been postulated that LiHoF4 is an example of 
an Ising-like system, and when the sample is placed in a 
magnetic field transverse to the Ising direction, it is an 
example of a quantum Ising system with magnetic dipole 
interaction being the ordering interaction. However, the 
quantum Ising model, Eq. Q has never been used to 
obtain physical results for the system. We derive the 
physical Ising model by a non-perturbative mapping to 
the Hilbert space spanned by the ground-state doublet. 
We then do a quantum Monte Carlo simulation to obtain 
the phase diagram, also incorporating the domain struc- 
ture in the process. This is a step beyond mean-field 
theory and the calculation is sufficiently accurate that 
uncertainty in the predicted phase diagram is now lim- 



ited by uncertainties in the crystal field parameters. As 
a spin-off, we are able to compute the phenomenological 
exchange interaction parameter that modifies the phase 
diagram considerably in the classical regime. The hy- 
perfine interaction poses considerable problems in com- 
paring the phase diagram obtained from the quantum 
Monte Carlo simulations to the experimental data. We 
have made the approximation that the effects of the hy- 
perfine interaction can be completely recovered through 
a renormalization of the magnetic field. 

LiHoF4 is a material in which the magnetic quantum 
and classical phase transitions can be controlled with 
great precision. Neutron scattering studies have been 
done on LiHoF4 to obtain the spin-wave excitations in 
the system. By randomly replacing the magnetic Ho 3+ 
with non-magnetic Y 3+ ions, spin-glass behaviour has 
been observed. We believe that the existence of an Ising 
model that faithfully reproduces the physics in both the 
classical and the quantum regimes, will facilitate the in- 
vestigations of the interesting properties of LiHoF4 con- 
siderably. 
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A. Appendix: The Crystal Field Hamiltonian 

In LiHoF4 , the Ho 3+ ions have an unfilled shell 4/ with 
10 electrons. The Hunds' rules dictate that the ground 
configuration of a single Ho 3+ ion should be 5 I$,{S = 
2, L = 6, J = 8). If there were no interactions with the 
neighboring ions, the ground state of a single ion will be 
17-fold degenerate. However, the Coulomb interactions 
with the neighboring ions gives rise to an electric field 
that lifts this degeneracy. In general, this electric field 
depends strongly on the spatial symmetry of the crystals. 
In the simplest scenario, each ion is regarded as a point 
charge, and the spatial overlap of the wave function of 
an ion with its neighboring ion is neglected. This is the 
point-charge model for calculating crystal fields. 

The derivation of the crystal field electrostatic poten- 
tial takes into account the lattice symmetry, and is most 
simply expressed in terms of the well-known spherical 
harmonics. We shall not discuss the details of the deriva- 
tion here, but instead refer the reader to the article by 
M. T. HutchingSf24 and references therein. If the elec- 
trostatic potential at a point r — (x, y, z) near the ion of 
interest is denoted by V(x,y, z), the crystal field hamil- 
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tonian can simply be written as 

Tie 

V c = -\e\Y i V{x i ,y i ,z i ). (52) 

i=l 

In Eq. (|52|) . n e denote the number of electrons in the 
unfilled shell (10 in case of LiHoF4), and (£i,yi 7 Zi) are 
the quantum position operators of the ith electron. 

Expressed in terms of the individual electron operators 
the crystal field hamiltonian still presents a formidable 
problem, but it is much simplified if we consider the 
outer shell electrons in an appropriate coupling scheme. 
In case of 4/ electrons, the ionic states are usually ex- 
pressed in terms of the quantum numbers \L, S, J,J Z ). 
If we restrict ourselves to a single value of J (neglecting 
configuration mixing), the hamiltonian in Eq. I|52l) can 
be expressed in terms of Operator Equivalents^^ This 
eliminates the need to use a Slater determinant consisting 
of single-electron wave functions to evaluate the matrix 
elements of Vc- 

The Operator Equivalents are operators built out of 
the J operators that act on the (2J + l)-dimensional 
space determined by the value of J. However, they act 
only on the angular part of the wave function of the cou- 
pled system, and the matrix elements of the radial part 
of the wave function are usually incorporated as fitting 
parameters. 

The number of operators needed to completely deter- 
mine the hamiltonian in Eq. (|52|) and the rules for de- 
riving them, depend on the symmetry of the crystal and 
the ground state configuration of the ion. These rules are 
clearly explained by K. W. H. Stevens^ Here we shall 
just list the operators that have non-zero matrix elements 
in the configuration 5 Ig of the Ho 3+ ion in LiHoF4. 

The operators are usually denoted by two indices, I 
and m, corresponding to the same indices on the spher- 
ical harmonics they are equivalent to. There is an addi- 
tional index, C(S), for non-zero m, corresponding to a 
symmetric(anti-symmetric) combination of the underly- 
ing spherical harmonics Y™ and Y j ~ m . 

In case of L1H0F4, the relevant crystal field operators 
are: 

2 ° = 3J 2 -J(J+1) 
Ol = 35J 4 -30J(J+l)J 2 + 25,/ 2 
-6J(J+1) + 3J 2 (J + 1) 2 

Oi(C) = \{Jl + Ji) 



O 6 = 231J 2 6 -315J(J+1)J 4 + 735J 4 

+ 105 J 2 ( J + l) 2 J 2 Z - 525 J( J + 1) J 2 + 294J 2 
— 5J 3 ( J + l) 3 + 40J 2 ( J + l) 2 - 60J( J + 1) 

Oi(C) = |(J+ + Jl)(lU|- J(J + l)-38)+/i.c. 

Ot{S) = Jl)(llJ 2 2 - J(J + 1)- 38) +h.c, 

(53) 

where J + = J x + iJ y and J_ = J x — iJ y . 

Using these operators, the crystal field hamiltonian Vc 
can be written as: 

V c = B 2 O 2 +B2O 4 + B° 6 O 6 +Bi(C)Ot(C) 

+Bt(C)Oi(C)+B*(S)Oi(S). (54) 

The radial matrix elements of the electrostatic poten- 
tial is extremely difficult to compute accurately even in 
a point-charge model. They are, therefore, incorporated 
within the constants Bf 1 , known as crystal field parame- 
ters(CFP). The CFPs are generally used as fitting pa- 
rameters. In LiHoF4, for example, they are used to 
fit the crystal field spectrum to observed spectroscopic 
( j a ^- a( io i 22 J 27 i 28 anc L £ susceptibility measurements^ 

In all our calculations, we use the CFPs proposed by 
R0nnow et al& Their values (in K) are listed below: 

B° 2 = -0.696 

B\ = 4.06 x 10~ 3 

B% = 4.64 x 10~ 6 

B\{C) = 0.0418 

B|(C) = 8.12 x 10~ 4 

B%(S) = 1.137 x 10~ 4 . (55) 

These values of CFPs were obtained 6 by fitting the 
results of RPA spin-wave dynamics calculations to ob- 
served neutron scattering data, as well as to the two 
lowest energy levels of the crystal field spectrum, as ob- 
served in spectroscopic measurements? 2 ^ However, there 
are no estimates of the accuracies to which these param- 
eters are known. There was an earlier attempt to de- 
termine the CFPs by fitting to susceptibility datajii but 
there were very large error bars. Another attempt was 
made to determine the CFPs by fitting to spectroscopic 
measurements^ but an incorrect symmetry(Z?2rf) of the 
crystal was used in the theoretical calculations, instead 
of the correct one, S4. 
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